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We show here how RenormaUzed Perturbation Theory (RPT) calculations applied to the quasi- 
linear growth of the large scale structure can be carried on in presence of primordial non-Gaussian 
(PNG) initial conditions. It is explicitly demonstrated that the series reordering scheme proposed in 
Bernardeau, Crocce and Scoccimarro (2008) is preserved for non-Gaussian initial conditions. This 
scheme applies to the power spectrum and higher order spectra and is based on a reorganization of 
the contributing terms into sum of products of multi-point propagators. In case of PNG new con- 
tributing terms appear, the importance of which is discussed in the context of current PNG models. 
The properties of the building blocks of such resummation schemes, the multi-point propagators, 
are then investigated. It is first remarked that their expressions are left unchanged at one-loop 
order irrespectively of statistical properties of the initial field. We furthermore show that the high- 
momemtum limit of each of these propagators can be explicitly computed even for arbitrary initial 
conditions. They are found to be damped by an exponential cutoff^ whose expression is directly 
related to the moment generating function of the one-dimensional displacement field. This extends 
what had been established for multi-point propagators for Gaussian initial conditions. Numerical 
forms of the cut-off are shown for the so-called local model of PNG. 

PACS numbers: 



I. INTRODUCTION 

With the advent of precision measurements of the large-scale structure properties of the universe (e.g. [1-3]), it 
became important to develop semi-analytic tools that allow to accurately compute the large scale statistical properties 
of the cosmic fields for any cosmological model. Indeed, it has been realized that the simple linear theory is too crude 
for the precision one wishes to attain. This is the case for instance for the precise determination of the power spectrum 
at the scales of the baryonic acoustic oscillations (see for instance [4-6]). 

The Standard Perturbation Theory approach (see [7]) is however not very efficient in producing well behaved next- 
to-leading order terms to the spectra or bispectra. It leads to perturbative series that have poor convergence properties 
(see [8] for more insights). Various resummation schemes have been proposed over the last few years to overcome 
those difficulties while aiming at providing tools for performing such calculations in a systematic way [8-13]. We will 
be interested here more particularly in the so called Renormalized Perturbation Theory (RPT) approach developed 
initially in [8] and [13] where rewritings of the resummation series are proposed. Whereas RPT as presented in [8] 
makes use of the one-point propagator only, in [13] it was introduced the concept of multi-point propagators that 
allow an alternative scheme for the computation of spectra and bispectra. 

Such developments however, have been so far limited to cases corresponding to Gaussian initial conditions (but 
see [14] for a scheme that numerically solves the truncated integro-differential equations governing the dynamics in 
presence of PNG) and it should be stressed that this is not only for convenience. In a diagrammatic representation 
of the dynamics, as presented for instance in [7] or in [8] the assumptions on the statistical properties of the initial 
field are indeed crucial. At the same time, new ideas to constrain the non-Gaussian properties of the initial metric 
fluctuations from observations of the large-scale structure of the universe have emerged (see for instance the recent 
reviews [15, 16] and [17]). They suggest that future redshift surveys could be our best opportunity to probe a departure 
from Gaussian initial conditions. This is in fact our main motivation for the investigations presented in this paper, 
where we extend the schemes and results flrst introduced in [13] to the case of primordial non-Gaussianities (PNG). 
The flndings presented here are rather general and, except for some numerical illustrations, are valid for any types of 
non- Gaussian initial conditions. 

The paper is organized as follows. In Sec. II we recall the equations of motion for gravitational instabilities of 
a cosmic fluid. In Sec. Ill we extend the expansion of power spectrum and bispectrum in terms of multi-point 
propagators to the case of arbitrary initial statistics. In Sec IV we compute explicitly the multi-point propagators in 
the large-A: limit for non Gaussian initial conditions. Lastly, in Sec. V we study the concrete example of primordial 
non Gaussianities of the local type. Our concluding remarks are presented in Sec. VI . 
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II. THE EQUATIONS OF MOTION 

We are interested here in the development of cosmological instabihties in a cosmological dust fluid. In general the 
dynamical evolution of such a fluid can be described with the Vlasov equation. As usual we restrict our investigations 
to the regime where multi-flow regions play a negligible role. In the one flow limit, the equations of motion then take 
the form of a set of three coupled equations relating the local density contrast, the peculiar velocity field and the 
gravitational potential (see [7]). 

At linear order these equations can easily be solved for an arbitrary background cosmology. One generically finds 
a growing solution and a decaying solution in the density contrast and peculiar velocity fields satisfying, 

^5iKT) + e(k,T) = 0, (1) 

1 d 0{k,T) / 1 d'H\9{k,T) 3^ ,,,,, , 



ndr H \ dr J H 2 

with (5(k, r) and ^(k, r) being respectively the Fourier transforms of the density contrast, (5(x, t) = /)(x)/p— 1, and 
of the pecuhar velocity divergence, = V • v. In Eqs. (1,2), Q,m{T) is the matter density and % = dlna/dr is the 
conformal expansion rate with a(T) the cosmological scale factor and r the conformal time. If one denotes D^{t) the 
growing mode solution of this system and /+ (t) its logarithmic derivative with respect to the expansion then, 

S{k, n) = D+(7?)(5o{k) , ^{k, n)/H = -f+{n)D+{rj)So{k) , (3) 

is the solution for the growing mode and 

S{k,rj) = D_{rj)So{k) , e{k,rj)/n = -f-{v)D_{rj)5o{k) , (4) 

for the decaying. We can note that in the Einstein-de Sitter limit we have /+ = 1 and /_ = —3/2, leading to = a 
and D_ = a~^/^. However, in the more generic case of arbitrary matter and dark-energy content and /_ do not 
admit an analytic solution [7]. 

Following [18], the equations of motion can be written in a compact form with the use of the two component 
quantity ^^(k, r), defined as 



MKr)^[5{k,T), -j-^0{k,T)), 



(5) 



where the index i = 1,2 selects the density or velocity components and which makes explicit use of the growing 
solution. It is then convenient to re-express the time-dependence in terms of the growing solution and in the following 
we will use the time variable rj defined as 

r7 = logD+(r7), (6) 

assuming the growing factor set to unity at initial time. Then the fully nonlinear equations in Fourier space read [7] 
(we henceforth use the convention that repeated Fourier arguments are integrated over). 



dr] 

where 



^^^(k,7y) + l^i,(7y)^,(k,77) =7^,-fc(k,ki,k2) *i(ki,7y) ^k{^2.v), (7) 



-1 

3 ^rn 3 1 

"2 fl 2 fl ^ 



(8) 



and the symmetrized vertex matrix ^y^jk describes the non linear interactions between diff"erent Fourier modes. Its 
components are given by 

7222(k,ki,k2) = (5D(k-ki-k2) --j-^ , 

7l2l(k,ki,k2) = ,5D(k-ki-k2) ^^^^^j^2'^\ (9) 
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with 7ij/c(k, ka, = 7i/cj(k, k5, k^), and 7 = otherwise, where denotes the Dirac delta distribution. The matrix 
^ijk is independent on time (and on the background evolution) and encodes all the non-linear couplings of the system. 
The formal integral solution to Eq. (7) is given by (see [8, 18, 19] for a detailed derivation) 



(10) 



where 0i(k) = ^i(k, = 0) denotes the initial conditions, set when the growth factor = 1 and where gij{r]) is the 
linear propagator^ that is the Green's function of the linearized version of Eq. (7) and describes the standard linear 
evolution of the density and velocity fields from their initial state. It is such that Qabiv^v') = for rj < r]^ due to 
causality, and Qahiv^v') ~^ ^ab as 77 — 77' ^ 0+. It naturally encapsulates the linear solutions growth and reads 





3 2 


D-iv) 
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f+iv) 



(11) 



for 77 > assuming the initial conditions (for a = 1) are set at a time when the universe is very close to Einstein-de 
Sitter. We note that growing mode initial conditions correspond to setting 0i(k) = 02 (k), in such a way that the 
second term in Eq. (11) does not intervene. The general expression (see [11]) of the propagator gij{ri^r]') is obtained 
from the property 



so that 

9ij{v,v') =9ik{r],^)g^j{r]',^). 
We can remark that for a Einstein-de Sitter background ^ij(ry, 77') = ^^^(77 — 77', 0). 



(12) 



(13) 



III. STATISTICS WITH NON-GAUSSIAN INITIAL CONDITIONS 



A. The F— expansion 



We are interested in the statistical properties of the density and velocity divergence fields, in particular in the 
construction of their power spectra and higher-order correlators. The (equal time) power spectra Pij are defined as 

(*,(k)*,(k')) = (5D(k + k')P,,(fc) , (14) 

and we want to reconstruct them from the expression of ^^(k) in terms of the initial density field 0j(k), the correlation 
properties of which are assumed to be known. The aim of this paper is precisely to explore the effects of dropping 
the assumption that the initial conditions are Gaussian distributed. 

A perturbative solution to Eq. (10) can be obtained by expanding the fields in terms of the initial ones, 

CO 

*,(k,r7) = ^*f)(k,,7), (15) 

n=l 

such that, 

= J d^qi...d^qn S^iq - qi...n) J^i^^j^..jjcii, • • • ^qn'^v) </>ji (Qi) • • • 0jn (^n) , (16) 

where we adopt the notation qi...n = Qi + • • • + Qn for vectors sums and where J^^^^ are fully symmetric functions 
of the wave- vectors that can be obtained recursively in terms of gij and 'jijk [7]. Note that these functions have a 
non-trivial time-dependence because they also include sub- leading terms in r]. Their fastest growing contribution is of 
course given by the well known {F^, Gn} kernels in PT (assuming growing mode initial conditions) whose dependence 
on time we recall here for an Einstein-de Sitter background, 

J^j^""^ = exp {nr]) {Fn(qi, Qn), Gn(qi, Qn)} , (17) 
for i = 1, 2 (density or velocity divergence fields respectively). 
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A formal expression for Pij{k) can be written using the expansion in Eq. (15) as, 

<5D(ki+k2m,(fci,r?) = ^ (*("^\ki,r?)*5"^)(k2,r?)). (18) 

ni ,n2 

For simphcity in what fohows we wih simply drop the component indices i and j but they are implicit. Then, for a 
given choice of indices rii and n2 one has to compute the ensemble average of rii + n2 factors ^(q*), following the 
field expansion in Eq. (16). The ensemble average of ^(qi) . . . ^(qni)^(qni+i) • • • ^(qni+n2) is given by the sum of 
product of cumulants of all set of subsets of {^(qi) . . . ^(q^ J^(qn^+i) . . . <p{q_n-^^-^n2)} ^^lat form a partition. Let us 
be more precise. Let us define Vn the set of partitions of a set of n variables ^. Its elements Si are lists of subsets, s^, 
and each element of Si is an index in the 1 . . . n range. As a result, 

(</.(qi)...</.(q„J^(q„,+l)...0(q„,+„j)= ^ n(n'^(q«)>c' (^9) 

where ( • • • )^ are the cumulants. Contrary to the case of Gaussian initial conditions, there exist cumulants that 
involve more than two variables. We can still assume though that there are no singleton although this hypothesis 
does not change the end of the calculation. 

The idea is now to sort the elements of Vni-\-n2- Each element Si of Vni-\-n2 defines four cardinal numbers; the 
numbers qi and q2 of points in subsets that are entirely within the first ni points or the last n2 points and the numbers 
Pi = rii — Qi and P2 = ^2 — Q2 of points within the first rii or last n2 that are in parts that are neither a subset of 
{1, . . . , ni} nor of {ni + 1, . . . , 712}. The resulting value of the partition contribution to the moment does not depend 
on which of the pi or p2 points is thus chosen because the functions are fully symmetric in their arguments. Up 
to symmetry factors (^j) {^l), it is then possible to assume that the pi points correspond to the first ones of ni and 
similarly for p2 . 

Let us now define 7^gi,pi,g2,P2 ^ subset of T^rii^n2^ ^-S- the set of partitions in T^yi^j^yi^ with fixed values of qi,Pi,q2 
and p2 that are formed from the union of an element of Vq^^ an element of Vq^ and an element of Vp^ p^. The latter 
is defined as the set of all partitions of {1, . . . ,pi + ^2} that do not contain subsets contained entirely within the first 
rii points or within the last n2. In Fig. 1 we give an explicit example of this partitioning scheme where an example 
of two different partitions having the same value are presented. 




FIG. 1: Example of partitions of ni + 712 points (here ni = 13 and 712 = 10). On each ensemble the elements have been 
numbered. For this choice of partitions there are 5 elements out of 13 that are in subsets that span both column and 4 out of 
10; therefore pi = 5, g'l = 8 and p2 = 4, ^2 = 6. The left hand side shows a generic partition of that kind. The right panel 
shows a partition which is identical to the left one up to a renumbering of the indices. These two contributing diagrams take 
the same values. When the order of the indices in each subset is preserved, the number of such contributing terms (e.g. the 
number of diagrams that lead to the same right hand side figure) is (^^^^^) (^^p^^^)- The summation over all possible partitions 
can then be restricted to those of the type of the left hand side with a weight given by the aforementioned symmetry factor. 
Such terms are defined by three new partitions that are elements respectively of Vq^ , Vq2 and Vp^ (^^^ text for definitions) . 



There is no simple expression for the number of such partitions. Its number is caUed the BeU number, Bn, after E.T. Beh, [20]. 
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As a result, the power spectrum in Eq. (18) reads 

Sb {^12) Pi jih^V) = / d^qi • • •d^qni+n2^D(ki-qi...ni)^D(k2-qni + l...ni+n2) X 

ni,n2 

^("^)(qi,...,q„J^("^)(qi+„„...,q„,+„J ^ H (H '^(*li))c ' (20) 



with 



where we have introduced the set of wave numbers {r^} that corresponds to reindexation of some of the q^, 

{ri, . . ..Vq^^q^} = {qi, . . . ,qg,,qn,+i, . . . ,qrii+q2}^ (22) 

and where Y^q^^p^^q^.p^ stands for =o Ep' =o ' Eg^^o Ep'^o' • The crucial property is then that the first two 
parenthesis that appear in the last line of the previous expression depend only on qi or q2 respectively. The sum over 
partitions can then be reorganized by summing for fixed values of pi and p2 first (see Fig. 1, right panel). That is, by 
doing 

n n-q . . / i \ 

e-^^"^ee(;;)-ee(';1^^^^^^ 

n>0 q=0 p=0 ^-^^ p>Oq>0 \ ^ / 

in each index 1 and 2 (after using that n = p -\- q) one can identify the function F^^^ defined as, 

r(f)(qi,...,qp,r7) = E ^ / d^qi . . . d^q, ^(f+«)(qi, . . . , q^, q^+i, . . . , qp+,; ,7) 

q=0 V 9 / ^ 

E n (u^i^p+i))]^ (23) 

(here i G [1,^]) which naturally extends to the case of arbitrary initial statistics the results for F^^^ studied in [9, 13] 
for Gaussian initial conditions. In such case q takes only even values and the two point initial spectrum determines 
the partition Vq, see for instance Eq. (15) in [9] and Eq. (21) in [13] for F^^^ and F^^^ respectively. 
After inserting Eq. (23) back into Eq. (20) we arrive at, 

fc(ki2)P(A:i,77) = X] / d^qi...d^qp, / d^q^.+i . . . d^q^.+p^fclki - qi...pi)fc(k2 - q^i+i.^pi+pj 

Pl,P2 

F(^^)(qi,...,q^,,77)F(^^)(q^,+i,...,q^,+^„7?) ^ JJ < 11 ^(^^)>c' (^4) 

^ Pl ,P2 

which is the final expression for the F-expansion of the power spectrum. It naturally extends the F-expansion obtained 
in [13] for Gaussian initial conditions to an arbitrary initial statistics. Note that contrary to that case the sum is not 
restricted to pi = P2- As a consequence there is in general no guarantee that all terms of this sum are positive. In 
addition, note that from Eq. (20) to Eq. (24) ^D(k— qi...p) could automatically be factorized out since ( Hien ^(%+*))c 
are all proportional to ^D(qp+i...p+q)- 

In Fig. 2 we show this expansion diagrammatically up to one-loop terms. Explicitly, these diagrams corresponds to 



P(k) = 





i> (1 

1) (2 






FIG. 2: F-expansion for the power spectrum up to one-loop. The shaded regions indicate the initial (or "primordial") statistics 
corresponding to each diagram, that is: one initial power spectra (2-points) in the left most diagram, one initial bispectrum in 
the center-left one, two power spectra in the center-right and one initial trispectrum in the right most figure. 



P{k) = [T^'\k)]^Po{k)^2T^'\k) yd3gF(2)(q,k-q)5o(q,k-q,-k) + 
+ 2 I d^g[F(2)(q,k-q)]2Po(k-q)Po(q) + 

+ J d^^id^g2r(2)(qi,k-qi)F(2)(q2,-k-q2)To(qi,k-qi,q2,-k-q2). (25) 

We stress that each term in the F-expansion above involves an infinite number of perturbative contributions in S 
as the F^^^ functions correspond to a full resummation of the propagator and the vertex in the language of [4, 10]. 

The resummation leading to the F-expansion for P{k) can obviously be extended for higher order correlators: higher 
order multi-point spectra can be obtained by gluing together F*^^^ functions multiplied by the proper cumulant. In the 
case of 3-point statistics (bispectrum) there are 3 "external" legs and the equivalent of Eq. (18) for the bispectrum 
P(ki,k2,k3) now runs over indices ni, n2 and 77,3. Following the same line of reasoning as above, the total set 
of partitions P^i+ns+na is broken up into three subsets Vq. contained fully with rii and one with "cross" elements 
^pi,P2,P3' bispectrum is then a sum over the product of elements of ^^^^^^^3 times F^^^^) F^^^^) pbs)^ 

as in Eq. (24). 

Figure 3 shows all the contributions in the multi-point propagator expansion of P(ki, k2, ka) up to one-loop diagrams. 
For PNG models satisfying the hierarchical scaling Bq ^ Pq ^ Tq ^ Pq and so on these are all the terms up to 0{Pq). 
These diagrams yield the following explicit expressions, 

B{h,k2,k3) = 2r(2)(ki,k2)r(i)(fci)r(i)(fc2)Po(fci)Po(fc2)+2perni.+ 

r(i)(fci)r(i)(fc2)r(i)(fc3) Bo(fci, fe, fcs) + 

8 J d3gr(2)(-q, q+ki) r(2)(-q-ki, q-k2) r(2)(k2-q, q) Po{q) Po(|ki+q|) Po(|k2-q|) + 

6 Po {ki )T'^'\ki)J d\ r(3) (ki , k2 - q, q) r(2) (k2 - q, q) Po (9) Po ( |k2 - q| ) + 5 perm. 

2Po(fci)r(i)(fci) r(2)(ki,k2)y"d3gr(2)(q,k2-q)Po(fc2,9,|k2-q|) + (fc2^fc3) +2perm. 

AT^^\k^)jd\ r(2) (q, k2 - q) r(2) (ki +q, k2 - q) Po ( |k2 - q| ) Bo (fci , |ki +q| ) + 2 perm. + 

3r(i)(fci)r(i)(fc2)Po(A:i) j d3qr(3)(ki,q,k2-q) 50(^2,*/, |k2-q|) + {ki ^ ^2) + 2perm. + 

r«(fci)r(i)(fc2)y"d39r(2)(q,k3-q)To(ki,k2,q,k3-q)+2perm., (26) 

where in the r.h.s. we implicitly assumed ka = — ki — k2 (and so the l.h.s only depends on the magnitudes of the 
wave- vectors). The first two contributions are "tree level". The next two correspond to the gravitationally induced 
non-Gaussianity, while the last four arise from PNG. 

Comparing with the standard perturbative approach up to one-loop given in [21] we see that his terms 5fi2, ^iii, 
^222 5 ^123 5 ^122 5 ^122 5 ^113 5 ^112 corrcspoud to thosc in Eq. (26) when the F^^^ functions are taken at their lowest or 
"bare" order in Eq. (23) (i.e. prior to any resummation). In turn, his ^fla (^113) correspond to the next-to-leading 
order term in the resummation of F^-*^^ in the first (second) line of Eq. (26). Lastly, his ^114 resums into F^^^ of the 
first line. 
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FIG. 3: F-expansion for the bispectrum up to one-loop. The first two contributions are proportional to P§ ("tree level"), the 
remaining ones to Pq , assuming that the bispectrum Bq and the trispectrum To scale respectively as Pq and P§. The meaning 
of shaded regions is described in the caption to Fig. 2 (see text for details). 



B. Multi-point propagators 



In the previous section we showed how the F^^^ functions in Eq. (23) serve as basic building blocks out of which 
one can construct series expressions for the polyspectra. In [13] these functions were further identified as n-point 
propagators since they can be obtained by functional differentiation of the final fields with respect to the initial ones. 
For PNG the F^^) functions are still equal to the ensemble average of the p-order functional derivative of ^ with 
respect to ^, that is 

pi (w£B^? ^ - ^s...- (1. • ■ ■ ^ ■ ") • (2') 

where we have re-introduced the doublet indices for clarity. This can be easily checked starting from Eq. (16) and 
getting to Eq. (23). Notice that, contrary to the standard FT kernels in Eq. (17), the multi-point propagators depend 
on the statistics of the initial fields. 



C. Multi-point propagators and correlation functions 



Within the framework developed so far, it is clear that a precise description of the large scale clustering can be 
achieved provided with a good understanding of the F^^^ functions. For Gaussian initial condition it has been observed 
that the F-functions can be expressed in terms of correlation functions between initial and final fields. For instance, 
the nonlinear propagator F*^^^ satisfies [9], 

(*i(k,7?)0,(k')) = (5D(k + k')r«(fc,r?)Po,;i(fc), (28) 

where Pq is the power spectrum of the initial density field and F^^^ is defined formally through functional differentiation, 
as in Eq. (27). An analogous expression is given in [13] for F^^^. These relations played a key role for the Gaussian 
case since they allowed the measurement of the fully nonlinear F^^^ and F^^^ function in N-body simulations [9, 13]. 

These simple expressions however, are no longer valid for non-Gaussian initial conditions. In this case the cross- 
correlation function between ^ and (j) can expanded as follows. 



(vl/(k, 7?) ,^(k')) = E / ■ ■ ■ ^D(k - qi...p) r(f)(qi, . . . , q^, ■ ■ ■ </'(qp)<^(q'))c. 



(29) 



This is actually a peculiar case of the general expression as obtained in Eq. (24) where the sum is now restricted to 
P2 = 1. This expression can then be written in terms of the spectra Po{k)^ bispectra 5o(ki,k2) and so on, of the 
initial fields times the corresponding F^^^ function from Eq. (23). It leads to, 

(vl/,(k,7?)(/.,(k')) = <5D(k-k') [r(/)(fc,r,)Po,y(fc) + 

d^qid^qz <5D(k-qi2) r|f^(qi, qz, v) So,imi(qi, q2) + ■ ■ • 1 (30) 
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The terms appearing in the previous equation are ordered in terms of importance for weakly non-Gaussian models 
such as the local [22-25] and equilateral [26] ones. Typically, the initial curvature bispectrum in these models satisfies 
the hierarchical^ scaling B^(k^ /c, k) ^ /nl P^{k), where P(f){k) is the curvature power spectrum. The second term in 
Eq. (30) is thus sub-dominant by a factor ~ J^lD^z). Alternatively, notice that Eq. (30) resembles the expansion for 
P(A:, 7^) in Eq. (25). This expansion have been studied using standard perturbation theory to one-loop and the leading 
order induced by primordial non Gaussianity (i.e. the second term in the previous equation) introduces corrections 
of ~ 1% — 2% for both local and equilateral models with \ fNL\ < 300 [27]. 

Similarly the computation of (^i(k, 77) (/)j(ki) (/)/(k2)) will make intervene the whole range of initial spectra. Its 
leading contribution is 

(^,(k, 7^) (/),(ki) (/)z(k2)> = ^D(k-ki2) [ Vf^{k, T]) 5o,m,7(kl, k2) 

+ 2 r,^^, (ki , k2 , ry) Po,mj {h ) Po,ni (0, ^^2) + • • • ] . (31) 

The relative importance of the terms in this series can be studied for particular models but in general, for weakly 
non-Gaussian initial conditions, the first two terms are dominant over the remaining series (not shown). Again, this 
expansion resembles the F-expansion for 5(/c, /ci, /C2) in Eq. (26). In turn, the bispectrum has been studied using 
one- loop perturbation theory by [21] for local and equilateral models, and additionally using N-body measurements 
in [28] for local models with /nl = ±100. In these cases, the next-to-leading terms become increasingly important 
for k ^ 0.1 Mpc depending on the value of /atl and redshift. 

Therefore, if we neglect sub- leading terms it is still possible to estimate the F- functions from correlators. This is 
well justified for T^^\ obtained from Eq. (30) as, 

^D(k - ki) r(i) {k, r,) = (*(k, r?) ,/.(ki)>/Fo(fc) • (32) 

Inserting this into Eq. (31) one recovers F^^^ as, 

A (^. v ^r(2)^i. v ^ (^(k,r;)0(ki)0(k2)) Po(ki, k2)(^(k, r;) 0(ki + k2)) 

2 Po (^1 ) Po (^2 ) 2 Po (^) Po (^1 ) Po (^2 ) 

As expected, the first term is the same as for Gaussian initial statistics. One can now insert this back into Eq. (30) 
and recover the correction introduced by the initial higher-order correlators such as Pq, 



wi, i^^rWrz. A (^(k,r;)0(ki)) ^ 1 /-^g,^ ^3,^ (^(k,r7)cA(ki)0(k2)) ^ ^ ^ ^ 



(34) 



In most commonly studied models of primordial non Gaussianity [22, 26], Pq ^ /nl- The expression above is thus 
neglecting terms of order /^^. In addition, as we discuss in the next section, the cross-correlations such as that in 
Eq. (34) will drop to zero towards high momenta. This is an important property to take into account when it comes 
to compute the momentum integration that cannot be recovered using the known tree-level analytical expression for 
these correlators. 

Expressions such as Eq. (33) and (34) give a concrete path to the actual measurements of the F functions that can 
certainly shed light into their full description. This is however beyond the scope of this paper and will be left for 
further work. 



IV. PROPAGATORS IN THE LARGE-fc LIMIT 



A. The a-method 



The multi-point propagators or F functions play an essential role in our formalism (as they did in [9, 13] for Gaussian 
initial conditions) and setting the ground for their adequate analytical description is thus one of the main goals of this 
paper. This task is however extremely difficult as it involves summing over infinite terms in Eq. (23) or adding up 



^ The extension of this scahng to the density spectra is however not straightforward because of the shape of the transfer function that 
relates the density and the potential at linear order. 
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an infinite number of diagrams, if we were to work with diagrammatic techniques as in [8, 13] (but see also [10, 11]). 
Nonetheless well established physical arguments were put forward in [9] to show that indeed, the leading order set of 
contributions can be resummed in the large-/c regime where the coupling of Fourier modes simplifies (see below). The 
low-A: regime, that allows perturbative calculations, could then be matched to the high-Zc asymptotic using simple 
physical arguments [9] . This resulted in a good agreement with propagator measurements in N-body simulations and 
allowed accurate descriptions of power spectrum [4] and bispectrum [13] at large scales. 

In what follows we are interested in extending these results to the case of non- Gaussian initial conditions, computing 
of the r functions (defined as the ensemble averages of functional derivatives of ipi with respect to the initial field 
values 0j) in the large- /c limit. As shown in [9, 13] this is a regime where the other intervening modes q in the 
momentum (or "loop") integrals are such that q <^ ki for all "external" k^. The intervening modes are also assumed 
to be in the linear regime with initial condition set in the growing mode. This later condition amounts to resum only 
those diagrams that are maximally connected with the initial conditions (that is, with all mode- mode interactions 
happening along the "principal path" as introduced in [9, 13]) what gives the leading contribution after resummation. 

Following these assumptions we then compute 9^^(k, 7^)/9(/)j(k') from Eq. (10) with the prescription that in its 
r.h.s. one replaces one factor by its linear solution ^/^/(t]', 0)0/. Furthermore 0/ is set is the growing mode, that is 
equal to So{q) ui with ui = (1, 1). From Eq. (10), one then has 



,(k') 



g,j{rj, 0) 5u{k - k') + 2 J\'D+{v') 9ik{v, v') lit (k, q, ki) m 5o(q) ^'^^^^^^l^,'^^ , (35) 



where the vertex function should be further taken in the q <^ ki limit (and recall that repeated Fourier arguments 
are assumed to be integrated over). As noticed in [9] this considerably simplifies its expression and we are left with 

^ 9«(,.0)*o(k -k') + io(k)/Vl'+W) V) . (36) 

where the intervening modes were all collected in a single quantity, 

a{k) = -iJd^ci^6o{q), (37) 

where So{q) is the initial density contrast. 

This possibility of collecting all the intervening modes in a single random variable (or more generally a finite number 
of variables) is at the heart of what we call here the a-method. This method was described and extensively used in 
[29] (and was already sketched in [30]). This method is much more powerful than standard diagram evaluations and 
countings. It eventually allows to compute the resulting propagator shapes either when the initial statistics is non 
trivial as in the present work and/or when the dependence of the propagators is complicated (as in [29]). 

Let us go back to this specific calculation. We note that for a given value of k, a(k) is a real random variable. It 
is Gaussian distributed if the initial conditions are Gaussian distributed but we will not make this hypothesis here. 
In the previous system we have introduced the i factor to make a(k) real and of positive variance. Note that in all 
cases, 

(a2(k)) = J d\ Po{q) ^ = fcVLp,., (38) 

where (Tdispi. is the one-component displacement dispersion in linear theory. 

We are then left with an equation evolution for d^i{k^a) / d(j)j{ri') that can be solved explicitly. If we indeed write 

^11^ = ^(k, r^) g,j{r, 0) fc(k - k') , (39) 

taking advantage of Eq. (12), we are left with the ODE 

— ^(k,r7)=ia(k)i?+(r)e(k,r), ^{k,rj = 0) = l, (40) 

which can be solved explicitly to give, 

= exp[ia{k){D+{rj) - lU^{r,,0)6^{k - k') . (41) 
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This is the first result of this section. We wih see later its implication for the form of the propagators. 

We can then proceed to the computation of higher order F functions from successive differentiations of Eq. (10) 
with respect to (j). The second order partial derivative induces two terms. 

d4>,{\.,)d<t>k{\^,) -y/^ 54^,^)7z„„(k,ki,k2) ^^^^^^^ ^^^^^^ 

^^J\'D,irj')Mv,v'Ut(^^^^^^^ (42) 
which can once again be simplified in the high- A: limit into, 

+ia(k) / drj^D^in^) 9uM) ^, .'LlV y (43) 
It can easily be checked that the solution of Eq. (43) is 



exp[ia(k)(D+(r7) - 1)] r dry' gu{v. v') itLi"^. ^iM) QmAv' . 0) QuM , 0), (44) 
Jo 

where the integral in the left hand part of this equation is the expression of 9^^i(k, a) / d(j)j(ki)d(j)k{k2) for a vanishing 
value of a that is its tree order expression. Indeed, inserting this expression in the second term of the right hand side 
of Eq. (43) and using the result of Eq. (41), we are left with 

52*,(k,r,) _ /"V ..... .'Us) 



/V 9ii{ri, n'hlt^iK ki, ka) exp[ia(ki + \^2){D+{ri') - l)\gmj{'n' , 0) QuM , 0) 



d(j)j{ki)d(j)kO<^2) 

ia(k) rdv'D+iv') ftz(??,r?')exp[ia(k)(I)+(V)-l)] /V gigiv' ,v") l^f^niKi^iM) gmjiv" ,0) gnkiv" ,0). (45) 

Jr]" Jo 

Then remarking that 

ia(k) /"%'D+(7?')ffa(»?,?7')exp[ia(k)(D+(?7') - 1)] gi^iv' = 

Jf]" 

5i,(7?,r/"){exp[ia(k)(Z)+(r,)-l)] -exp[ia(k)(Z)+(7?')-l)]}, (46) 

and that a(k) is a linear function of its argument, the two terms of the right hand side of Eq. (45) simplify and 
recombine into the expression of Eq. (44) . 

This result can naturally be extended to any order^ as it can be established recursively by successive use of the 
relation (46). More precisely, we have 



-exp[i.(k)(i.,-l)] 



(47) 

a(k)=0 



dcj^j,{k,)...dcj>j^{kp) ' ^ - a(/),-,(ki)...a0,-^(kp) 

where k = X] kj . After taking the ensemble average of this expression we finally obtain, 

rW(ki, . . . , kp, 7?) = (exp [ia (k) (D+(,?) - 1)] ) r[±{ku . . . , k^, ,?) , (48) 

where we have used that the derivatives on the r.h.s of Eq. (47) are evaluated at = (i.e. at tree level) and 
thus are independent of initial random fields. This is a remarkable result that extends in a concrete and simple way 
the findings of [9, 13] for Gaussian initial conditions. Hence, the statistical properties of "random variable" a(k) 
determine a transition function that relates the fully nonlinear multi-point propagator with their tree level expression 
(which, up to sub-leading terms, are nothing else than the standard PT kernels in Eq. 17), in the high-A: regime. In 
what follows, we will use this result to explicitly compute the F- functions in this regime. 



^ An alternative way of obtaining this result is to remark that in presence of low q perturbing modes Eq. (7) can be rewritten - for 
an Einstein-de Sitter background - as a^^i(k, a) + r2^j^j(k, a) — iQ;(k)a^^(k, a) = 7^jfc(k, ki, k2) ^j(ki,a) ^fc(k2,a). The linear 
progator is (k, ai , 02) = exp(iQ;(k)(ai — a2)) Qij (ai / a2) • The quantities we derive in this paper are the tree order F-functions 
of this theory, whose a dependence can be easily computed from the fact that ^ij(ki + k2, ai, a2)^A;z(ki, a2, a3)^rnn(k2, a2, as) = 
exp(ia(ki +k2)(ai - a^))gij{ai/a2)gkl{0'2/ci^)gmn{ci2/ci^)- 
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B. The transition function 

Irrespectively of their order, the F-functions are all obtained, in the large-A: regime, from the computation of the 
ensemble average in 

/(fc) = (exp[ia(k)(D+-l)]), (49) 

that will depend only on k due to homogeneity and isotropy. From the definition of a(k) in Eq. (37) we observe that 
the function f{k) can in fact be easily related to the statistical properties of the primordial field: it is nothing but 
the moment generating function M (see Eq. (133) of [7]) at point t = k(D+ — 1) of the single point one component 
displacement field d. 



d = J d3q^,5o(q). 



(50) 



As a consequence of Eq. (135) in [7], log f{k) is nothing but the cumulant generating function of that same variable 
at the same point defined through, 

iog/(^)=f:^^^^(^+-i)^- (51) 

p=2 P' 

For Gaussian initial conditions only p = 2 is non-zero and one recovers the well known result f{k) = 
exp(— /c^cr^-gpj (I^+ — 1)^/2) [9, 13]. Note also that only even values of p contribute to this sum because the cu- 
mulants cannot depend on the direction of the wave-vector k. This is at variance with the result derived in [31] 
following the prescription for propagator resummation put forward in [10] for Gaussian initial statistics. It is however 
a very important result because it shows that, for any given model of primordial non-Gaussianity, the first "non- 
Gaussian" correction to the propagator decay (in the high-/c limit) is given by the four point connected function. For 
weakly non-Gaussian initial conditions this will represent a minor contribution. 

V. PREDICTIONS FOR THE LOCAL MODEL OF PNG 

We now turn into the evaluation of the transition function f{k) defined in Eq. (51) for a specific model of primordial 
non-Gaussianity. We will focus in the local model [22-25], which is perhaps the most studied model of primordial 
non Gaussianity within the context of large scale structure (see for instance the recent reviews [15, 16] and references 
therein). This model is build upon the nonlinear relation [22] 

$(X) = (/.(X) + fNL{^\^) - {^H^m + QNL^H^) , (52) 

between the Gaussian field and Bardeen's gauge- invariant potential ^, where up to cubic terms are considered. The 
local model corresponds to /nl and g^L constants independent of space. The connection to the matter overdensity 
is through Poisson's equation ^k(^) = ^(^,^)^k, with 

«(M) = |«^^, (53) 

where T{k) is the matter transfer function, and D{z) the growth factor. In what follows we will assume the same 
cosmological model as in [28] (and the CMB convention for f^L and gNh)- This can be summarized as h = 0.7, 
Vim = 0.279, = 0.0462, = 0.96 and erg c± 0.81 (see Sec. 3 in [28] for more detail) 
We start by writing the first terms of the sum defining /(/c), 

log/(fc) = i^^^(D+ - If + ^^^^(Z)+ - 1)^ + . . . , (54) 

with the displacement field d given by Eq. (50). As noticed before, the first contribution reduces to ((d • k)^)c = 
~^disp^^ and (Jdisp = 6.01 Mpc (at 2: = 0) for our cosmology. For Gaussian initial conditions this is the only 
non-zero contribution. For non-Gaussian initial conditions we need to evaluate, in addition 

((d-k)^), = /dVfrQ2frg3/dV'^'^'^'^(^o(qi)^o(q2)<5o(q3)<5o(q4))c 

J J J J Q.I Q.2 % Q-A 

f^s f^s ^3 k qi k q2 k qs k .q4 

= d qi d q2 d qs — 2 2 2 2" ^o(qi, q2, Qs, q4) , (55) 

J J J Qi Q2 % ^4 
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k (h Mpc-i) 

FIG. 4: The non Gaussian transition function f(k) in the large- /c limit in ratio to its expression for Gaussian statistics (i.e. 
setting fNL^QNL^O). For currently accepted values of non Gaussian parameters (i.e. close to those shown) the correction 
to the damping is remarkably sub-dominant, even more so at higher redshift (e.g. bottom panel). 



with q4 = — qi23 = — (qi +^2+^3) in the last equahty. The matter trispectrum Tq is given by 

To(ki,k2,k3,k4) = M{ki,z)M{k2,z)M{ks,z)M{h,z) T^{ki,k2,ksM) . (56) 

where is the curvature trispectrum. For local non-Gaussianity there are two distinct contributions, 

T^(ki,k2,k3,k4) = 4/|,^P^(fei)P$(fe2)[P$(^i3)+P$(^i4)]+5perm. 

-^6gNLP^{ki)P^{k2)P^{ks) + 3 perm. . (57) 

Hence, in what follows we will distinguish the /^^ component from the qnl one. In appendix A we include a 
detailed account of the evaluation of the integral in Eq. (55). The final result for our cosmology is ((d • k)^)c = 
fNLk^^tjj^L + d^Lk^^XoNL w^^^ ^^Jnl = 0.05537 /i"^Mpc and 54,^^^ = 0.02359 /i"^ Mpc (at z = 0). After putting 
together this result and the one for the variance into Eq. (54), and normalizing the growth to 1^(2; = 0) = 1, we obtain 
the prediction for the transition function in the local model at arbitrary redshift. 

Figure 4 shows the ratio of the f{k) function so obtained, assuming Jnl = 10^ and/or Qnl = 10^, to its value for 
Gaussian initial conditions (i.e. with the kurtosis set to zero). We see that the modification to the Gaussian case 
for f{k) is minor for this choice of the PNG parameters, with up to 2% weaker decay at z = (top panel) and 1% 
at z = 0.5 (bottom panel). In turn, Fig. 5 shows log/(/c) as a function of k'^ for Gaussian and non-Gaussian initial 
conditions, with /nl^^x 10^ and, separately, Qnl = 10^. As expected, the correction becomes relevant for values 
of Jnl that make the (jp' correction the same order of the Gaussian component. 

Provided with the damping function f{k) we are in place to explicitly compute the F -expansion for the bispectrum, 
as given in Eq. (26), but assuming the multi-point propagators in their large-/c limit (a well justified approximation, 
at least for Gaussian initial conditions [13]). Figure 6 shows the terms contributing to this expansion assuming a 
local model for primordial non-Gaussianity with /nl = 300 and Qnl = 10^, and equilateral configurations. In solid 
blue line we show the equivalent to the tree- level bispectrum induced by gravity, i.e. the first term in Eq. (26). The 
contribution from the primordial bispectrum is depicted by solid red line, it dominates the total signal at very large 
scales but it is exponentially suppressed at high-/c by F*^^^. At one-loop there are two terms also present for Gaussian 
initial conditions, see discussion after Eq. (26), which are 0(Pq). These are shown by the dashed red line. The 
purely PNG induced ones at this order are given by the solid green line, they compromise terms 0{PqBq) ~ /nlPq 
and 0(To) ^ (/nl^ 5'nl)Po • dashed line we show the result of standard perturbation theory at one-loop from [28]. 
Contrary to that case, here the number of terms is reduced (due to the resummation of several contributions). In 
addition, each term dominates over a narrow range of scales. 
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FIG. 6: The F expansion for the equilateral bispectrum in the local model of PNG, up to one-loop. The "gravitational" tree 
level bispectrum is shaped by F^^-* while the PNG one arises from the primordial bispectrum and is exponentially suppressed 
to high-A: by F^^\ At one-loop there are terms Pi (grav.), ^ Pq x Bq (PNG, /nl) and To (PNG, /nl and ^nl). 



VI. CONCLUSIONS 



The impact of primordial non-Gaussian (PNG) initial conditions on the statistical properties of the cosmic density 
field is a priori not easy to infer. This is because most of the diagrammatic expansions that have been lately used to 
study the large scale clustering are rooted in the assumption of primordial Gaussian fluctuations, that considerably 
simplifies the structure of multi-loop diagrams. We show here explicitly that PNG can indeed be accounted for in the 
approaches developed in previous papers along the so-called RPT method [8] . In particular the F-expansion approach 
of [13] is preserved for such initial conditions. 

Within this context we generalize the definition of multi-point propagator to the case of arbitrary initial statistics. 
We find that they can still be regarded as the basic building blocks out of which one can construct n-point spectra. The 
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F-expansions so obtained correspond to a resummation of infinite sub-sets of diagrams in the approach of [8, 10, 12], 
concretely those corresponding to propagator and vertex renormahzation, in hne with the results in [13] for initial 
Gaussian fields. For concreteness we give explicit expressions for the series expansion up to one-loop of the power 
spectrum and the bispectrum. In this way, the leading corrections to the power-spectrum are then those that make 
intervene the primordial bispectrum in the contraction of the F^^^ — r^^^ product. 

We then turned into the description of the F*^^^ functions themselves (a.k.a. the multi-point propagators). We show 
how the simple one-to-one correspondence between n-point propagators and n-point correlators (among initial and 
final fields) found for Gaussian initial conditions can be extended to the case of PNG. In the latter case, however, we 
find that this relations involves a infinite hierarchy among F^^^ and correlators. Nonetheless for weakly non- Gaussian 
models, as those allowed by current data, well justified approximations can be put forward to close the hierarchy 
permitting the actual measurement of F^'^^ in N-body simulations. In this way one can extend the methodology 
carried out for Gaussian initial statistics [9, 13] to shed light in the full description of the multi-point propagators. 

Finally, the large- /c damping behavior of the multi-point propagators is derived for arbitrary initial conditions. In 
this regime, we find that the multi-point propagators are proportional to their tree- level expressions. Moreover we 
explicitly show that the rate of decay is the same irrespectively of their order, and that is given by the cumulant 
generating function of the displacement field. Remarkably, the first corrections to the "Gaussian" damping in this 
regime is due to the initial trispectrum, which for small departure from Gaussian initial conditions constitutes a 
sub- dominant contribution. 

We leave for future work a more detailed analysis involving quantitative predictions for specific models of PNG. 
But we hope the results already presented here might serve as a basis to meet the accuracy requirements demanded 
by the analysis of future datasets tailored at deciphering the primordial statistics and the physics of inflation. 
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Appendix A: Evaluation of the kurtosis cumulant 



In this appendix we discuss the evaluation of the 4-order cumulant defining the first non- Gaussian contribution to 
the transition function f{k) in Eq. (51). As discussed in Sec. V, and particularly in Eqs. (56,57), we must evaluate 
two separate contributions to the trispectrum defining the kurtosis cumulant. We have, in the first place 

((d . k^),,,, = U^q, ^^^^ r,„.(qi,q2,q3,q4) , (Al) 

J J J HI H2 % Q4 

with q4 = — qi23 and 

Tf^^ (qi, q2, qs, q4) = 48 /^^ M(gi) M{q2) M{qs) Miq^) P^{q^) P^{qs) P^{qi2) , (A2) 

giving the /nl term, and a similar expression with the trispectrum replaced by 

f^^ Jqi, q2, qs, q4) = 24gNL M{q^) M{q2) M{q^) Miq^) P^{q^) P^{q2) PM . (A3) 

corresponding to the qnl contribution. Notice that in Eqs. (A2,A3) we took advantage of the fact that we are 
integrating over all orientations so all permutations give the same result. Since both integrals depend only on the 
magnitude A:, we can integrate over the orientation of k. Hence, for the Jnl contribution we have 

((d-k)4),/„, = Jd^q^Jd%Jd%ff^A<li,ci2,cis,ci,)l- Jdn, (A4) 

and we are then allowed to put qi along the 2:- axis and set the azimuthal angle of q2 equal to zero, obtaining the 
6-dimensional integral 



((d • k)^)c,/^^ = 27rJdqiqlJdq2qlJdcos02jd^q3Tf^^{qi,q2,c[3^Ci4)Jdnk 



k • qi k • q2 k • q3 k • q4 

2 2 2 2 ' 
QI Qi Qi 



2^|dgi| dg2 |dcos^2 Jd^Qs ^^'^ g^f ^''^ Jd^k (k • qi)(k • q2)(k • q3)(k • q4). 



QsQi 

Introducing the vector p = qi2 = qi + q2 we can rewrite the integral as 

I 



{(d.k)^)e,/^, = 2n[ dqJ dq2r^'^dp^[ d'qs ^f-^^^^^ff^^^^^) f (k • qi)(k • q2)(k • q3)(k • q4), 

J J Jlqi-qol Q^Q'^J Q3Q4 J 



and recall that q4 = — qi23 = — qi2 ~ q3 = ~P ~ q3 5 so that the magnitude q^ only depends on the magnitudes of p 
and qs and on the angle between them. Then, we can transform the integrand in the basis where p = qi2 is along 
the 2;-axis, 

. P [a ^[a a Tf^^{qi,q2,q^,qA,p) 



((d • k)^)e,/^^ =2^/ dqJ dq2f ' ' dp^f dqsqjf dcos^s 

J J J\qi-q2\ Q1Q2J J 



qIqI 



J d4>3 J dQk (k • qi)(k • q2)(k • q3)(k • q4), (A5) 
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with Os and ^3 define the orientation of qs in the new basis and Os is now the angle between qs and p. Given the 
one-to-one correspondence between 63 and ^4, as before we can change variable of integration to get 



dqi dq2 dp^UqJ Aq, 

J J\qi-q2\ ^^(I'iJ J\p-qi\ 

/r rqi-tq2 r rP-tqs 

Mil M2I dp dqs dg'4 
J Aqi-q2\ J Av-q3\ 
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PQ3 qi 



dQk (k • qi)(k . q2)(k • q3)(k • q4) (A6) 

qiq2q3qA 

X Fgeom {qi , ^2 , gs , ^4 , P) , ( A7) 

with the geometric factor given by, 

-Fgeom (^1 •> q2,q3,qi,p) = J '^^^ J ^^fc (k-qi)(k-q2)(k-q3)(k-q4) 
= ^M-Qif-7{ql + ql)p'+6p'] 

M - & - 3(9? + iDp' + 2p1 
+ 1^ [3(9? - 4f - 5(9? + ql)p' + 2/] (A8) 

being a function of the five variables. The g^L integral follows from the same considerations, with T^^^(gi, qs, ^4) 
replacing Tj^^(gi, ^2, qs^q^^p)- An explicit evaluation of these integrals for the cosmology assumed in this paper (see 
Sec. V) yields 

((d • k)^)e = f^L 9-4 X 10-^ + QNL 3-1 X 10"^ (A9) 
assuming Tq is evaluated at 2; = and k is in units of /iMpc""^. 



